// STL
#include <stdexcept>                           // std::runtime_error()
#include <string>                              // std::string, std::to_string()
#include <tuple>                               // std::tuple<>, std::make_tuple(), std::tie()
#include <vector>                              // std::vector<>

// Boost
#include <boost/math/distributions/normal.hpp> // boost::math::normal

// jScience
#include "jrandnum.hpp"                        // rand_num_normal_Mersenne_twister()

// functions++
#include "functionsxx/regression.hpp"          // functionsxx::regression::Deming()

// stats++
#include "statsxx/postprocess/ROC.hpp"         // ROC_pt
#include "statsxx/statistics.hpp"              // statistics::mean(), ::stddev()


//
// DESC: Estimate (fit) an ROC curve to a binormal model.
//
// INPUT:
//
//     std::vector<ROC_pt> ROC         :
//     unsigned int        n           : (default of 100). Note the actual number of points stored in the ROC curve will actually be two less, because the end points are discarded.
//
// OUTPUT:
//
//     bool                success     :
//     std::string         msg         :
//     double              a           : binormal coefficient
//     double              b           : binormal coefficient
//     std::vector<ROC_pt> binormalROC : binormal curve (evaluated at fitting points)
//
// NOTE: Derivation:
//
//     Assume that the responses on positive (P) and negative (N) cases both follow normal distributions with means mu and variances sigma^2 ---
//
//     ... i.e., N(mu_p, sigma_p^2) and N(mu_n, sigma_n^2), respectively.
//
//     For an arbitrary point mu:
//
//          TPR = P(x_p > mu) = 1 - Phi((mu - mu_p)/sigma_p)
//          FPR = P(x_n > mu) = 1 - Phi((mu - mu_n)/sigma_n)
//
//     ... where Phi is the distribution function of a standard normal variable.
//
//     If follows, for example, that:
//
//          mu = mu_n + sigma_n*Phi^-1(1 - FPR)
//
//     Therefore:
//
//          TPR = Phi(a + b*Phi^-1(FPR))
//
//     ... (which does not depend on mu) where:
//
//          a = (mu_p - mu_n)/sigma_p
//          b = sigma_n/sigma_p
//
//     a and b can therefore be found by the linear curve defined by:
//
//          Phi^-1(TPR) = a + b*Phi^-1(FPR)
//
// NOTE: Bootstrapping is used to calculate Phi^-1(TPR) and Phi^-1(FPR) and associated errors.
//
// Notes:
//     - I ?think? that it is not possible to extract the *explicit* normal distributions, only the *latent* ones
//
//
// TODO: NOTE: This subroutines assumes that there is (independent) error in the x and y directions (by threshold averaging) --- both through the use of Deming regression, and the assignment back of points and ROC scores.
//
// TODO: NOTE: ... Should probably do checking of the input, or take more user input to figure this out.
//
inline std::tuple<
                  bool,
                  std::string,
                  double,
                  double
                  > fit_binormal(
                                 std::vector<ROC_pt> ROC,
                                 std::vector<ROC_pt> stddevROC,
                                 // -----
                                 const int           nROCs
                                 )
{
    bool                success;
    std::string         msg;
    double              a;
    double              b;

    //=========================================================
    // INITIALIZE
    //=========================================================

    boost::math::normal dist;

    // NOTE: The quantile function blows up in magnitude at 0 or 1.
    //
    // NOTE: ... Check that two asymptotic points [(0,0) and (1,1)] exist, and remove them.

    if( (ROC.front().FPR == 0.) && (ROC.front().TPR == 0.) )
    {
        ROC.erase(ROC.begin());
        stddevROC.erase(stddevROC.begin());
    }
    else
    {
        throw std::runtime_error("error in fit_binormal(): first point in ROC curve is not (0,0)");
    }

    if( (ROC.back().FPR == 1.) && (ROC.back().TPR == 1.) )
    {
        ROC.pop_back();
        stddevROC.pop_back();
    }
    else
    {
        throw std::runtime_error("error in fit_binormal(): last point in ROC curve is not (1,1)");
    }

    for( auto i = 0; i < ROC.size(); ++i )
    {
        if( (ROC[i].FPR == 0.) || (ROC[i].FPR == 1.) || (ROC[i].TPR == 0.) || (ROC[i].TPR == 1.) )
        {
            throw std::runtime_error("error in fit_binormal(): points on the ROC curve lie on (x,0/1) or (0/1,y). try threshold averaging with less points.");
        }
    }

    //=========================================================
    // FIT THE SHAPE OF THE ROC CURVE TO A BINORMAL MODEL
    //=========================================================

    //---------------------------------------------------------
    // COMPUTE Phi^-1(TPR) AND Phi^-1(FPR) (AND ERRORS)
    //---------------------------------------------------------
    //
    // NOTE: Bootstrap is used for this, in order to get the errors.

    std::vector<double> Phiinv_FPR;
    std::vector<double> Phiinv_FPR_stddev;
    std::vector<double> Phiinv_TPR;
    std::vector<double> Phiinv_TPR_stddev;

    for( auto i = 0; i < ROC.size(); ++i )
    {
        std::vector<double> _Phiinv_FPR;
        std::vector<double> _Phiinv_TPR;

        int n = 0;

        while( n < nROCs )
        {
            double _FPR = rand_num_normal_Mersenne_twister(
                                                           ROC[i].FPR,
                                                           stddevROC[i].FPR
                                                           );

            double _TPR = rand_num_normal_Mersenne_twister(
                                                           ROC[i].TPR,
                                                           stddevROC[i].TPR
                                                           );

            if( (_FPR > 0.) && (_FPR < 1.) && (_TPR > 0.) && (_TPR < 1.) )
            {
                _Phiinv_FPR.push_back( boost::math::quantile( dist, _FPR ) );
                _Phiinv_TPR.push_back( boost::math::quantile( dist, _TPR ) );

                ++n;
            }
        }

        Phiinv_FPR.push_back( statistics::mean( _Phiinv_FPR ) );
        Phiinv_FPR_stddev.push_back( statistics::stddev( _Phiinv_FPR ) );

        Phiinv_TPR.push_back( statistics::mean( _Phiinv_TPR ) );
        Phiinv_TPR_stddev.push_back( statistics::stddev( _Phiinv_TPR ) );
    }

    //---------------------------------------------------------
    // FIT: Phi^-1(TPF) = a + b*Phi^-1(FPF)
    //---------------------------------------------------------

    // TODO: NOTE: Deming regression assumes that there is (independent) error in the x and y directions.

    std::vector<double> Phiinv_FPR_star;
    std::vector<double> Phiinv_TPR_star;
    std::tie(
             a,
             b,
             Phiinv_FPR_star,
             Phiinv_TPR_star
             ) = functionsxx::regression::Deming(
                                                 Phiinv_FPR,
                                                 Phiinv_FPR_stddev,
                                                 Phiinv_TPR,
                                                 Phiinv_TPR_stddev
                                                 );

    // NOTE: By definition, [a = |mu_s - mu_n|/sigma_s] or [b = sigma_n/sigma_s] cannot be negative.
    if( (a < 0.) || (b < 0.) )
    {
        throw std::runtime_error( ("error in fit_binormal(): a or b is negative. a = " + std::to_string(a) + "; b = " + std::to_string(b)) );

//        std::string msg = "error in fit_binormal(): a or b is negative. a = " + std::to_string(a) + "; b = " + std::to_string(b);
//
//        return std::make_tuple(
//                               false,
//                               msg,
//                               a,
//                               b,
//                               std::vector<ROC_pt>()
//                               );
    }

    // TODO: NOTE: Maybe here add a warning or something is the binormal curve is improper.

    //=========================================================
    // RETURN
    //=========================================================

    return std::make_tuple(
                           true,
                           std::string(),
                           a,
                           b
                           );
}
